# Analysis of Streamflow data
# Created as part of Contract GSA-IS17JHQ0177
# Author Carl James Schwarz
# StatMathComp Consulting by Schwarz, Inc.
# cschwarz.statsfuca@gmail.com
# Change Log
# 2017-01-30 CJS First Edition
options(width=200) # how wide to print output
# load the functions used in the analysis
# This will also load the relevant libraries -- see the loadlibraries file in the Functions folder for details
dummy <- lapply(file.path("Functions",list.files(path=file.path("Functions"), pattern = "(?i)[.]r$", recursive = TRUE)),
function(x){
cat("Loading ",x,"\n");
source(x)})
## Loading Functions/compare.annual.stat.R
## Loading Functions/compare.frequency.with.hec.R
## Loading Functions/compare.longterm.stat.R
## Loading Functions/compare.percentile.longterm.stat.R
## Loading Functions/compute.Q.percentile.longterm.R
## Loading Functions/compute.Q.stat.annual.R
## Loading Functions/compute.Q.stat.longterm.R
## Loading Functions/compute.volume.frequency.analysis.R
## Loading Functions/load_packages.r
## Loading required package: MASS
## Loading required package: survival
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Define variables that provide options for the analysis
Station.Code <- '08HA011'
Station.Name <- 'Upper Cowichan River'
Station.Area <- 826 # square km's
start.year <- 1965 # when do you want the analysis to start at.
end.year <- 2012 # what is last year with complete data
# Get the data
flow <- read.csv(file.path("08HA011","08HA011_STREAMFLOW_SUMMARY.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
# Create a date variable
flow$Date <- as.Date(paste(flow$Year,flow$Month, flow$Day, sep="-"), "%Y-%m-%d")
# basic structure of flow
str(flow)
## 'data.frame': 42734 obs. of 7 variables:
## $ Day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Month: int 1 1 1 1 1 1 1 1 1 1 ...
## $ Year : int 1900 1900 1900 1900 1900 1900 1900 1900 1900 1900 ...
## $ DOY : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Q : num NA NA NA NA NA NA NA NA NA NA ...
## $ Hour : chr "24:00" "24:00" "24:00" "24:00" ...
## $ Date : Date, format: "1900-01-01" "1900-01-02" "1900-01-03" "1900-01-04" ...
# Some preliminary screening
# Are there any illegal dates that could not be created?
cat("Number of illegal date values ", sum(is.na(flow$Date)), "\n")
## Number of illegal date values 0
flow[ is.na(flow$Date),] # sorry Microsoft, but feb 29 1900 is NOT a leap year
## [1] Day Month Year DOY Q Hour Date
## <0 rows> (or 0-length row.names)
dim(flow) # size before removing illegal dates
## [1] 42734 7
flow <- flow[ !is.na(flow$Date),]
dim(flow) # size after removing illegal dates
## [1] 42734 7
# Table of simple statistics by year to help check for outliers, etc
flow.sum <- ddply(flow[ flow$Year >= start.year & flow$Year <=end.year,], "Year", plyr::summarize,
n.days = length(Year),
n.Q = sum (!is.na(Q)),
n.miss.Q = sum ( is.na(Q)),
min.Q = min (Q, na.rm=TRUE),
max.Q = max (Q, na.rm=TRUE),
mean.Q = mean(Q,na.rm=TRUE),
sd.Q = sd (Q,na.rm=TRUE))
flow.sum
## Year n.days n.Q n.miss.Q min.Q max.Q mean.Q sd.Q
## 1 1965 365 365 0 3.54 196.0 47.72625 44.55059
## 2 1966 365 365 0 5.80 362.0 60.32545 65.58664
## 3 1967 365 365 0 4.11 213.0 61.70992 53.65439
## 4 1968 366 366 0 5.10 450.0 69.52503 71.15098
## 5 1969 365 365 0 4.98 150.0 47.50655 35.52722
## 6 1970 365 365 0 4.25 153.0 36.84767 34.33385
## 7 1971 365 365 0 6.00 187.0 55.62266 39.62329
## 8 1972 366 366 0 5.83 425.0 53.93230 59.41003
## 9 1973 365 365 0 3.11 337.0 53.33545 64.81365
## 10 1974 365 365 0 3.28 343.0 69.46066 64.42318
## 11 1975 365 365 0 6.43 345.0 66.59186 65.38940
## 12 1976 366 366 0 4.62 165.0 47.77153 37.67202
## 13 1977 365 365 0 4.16 233.0 47.69573 47.43126
## 14 1978 365 365 0 4.64 101.0 30.25252 21.66569
## 15 1979 365 365 0 3.84 336.0 43.46542 56.16374
## 16 1980 366 366 0 4.42 385.0 58.65967 62.00653
## 17 1981 365 365 0 4.66 244.0 60.31227 54.73371
## 18 1982 365 365 0 5.32 217.0 49.76348 46.18497
## 19 1983 365 365 0 5.71 303.0 61.29674 69.73937
## 20 1984 366 366 0 4.12 197.0 58.12112 41.60186
## 21 1985 365 365 0 3.36 86.4 27.75455 20.36948
## 22 1986 365 365 0 3.45 400.0 52.85359 63.78841
## 23 1987 365 365 0 3.60 198.0 47.36142 50.49256
## 24 1988 366 366 0 4.14 139.0 46.11101 37.02558
## 25 1989 365 365 0 4.08 176.0 40.17173 35.15942
## 26 1990 365 365 0 4.89 326.0 59.71548 63.55819
## 27 1991 365 365 0 4.61 248.0 45.93523 47.93800
## 28 1992 366 366 0 3.74 324.0 42.74970 56.93322
## 29 1993 365 365 0 4.89 191.0 33.46789 33.87104
## 30 1994 365 365 0 4.80 257.0 47.67293 52.08854
## 31 1995 365 365 0 4.67 286.0 71.34422 77.38259
## 32 1996 366 366 0 4.63 254.0 50.60377 45.75725
## 33 1997 365 365 0 7.45 359.0 79.12233 55.64481
## 34 1998 365 365 0 3.96 325.0 57.35233 71.50327
## 35 1999 365 365 0 6.01 282.0 71.80063 57.92698
## 36 2000 366 366 0 5.23 88.7 36.60932 22.77473
## 37 2001 365 365 0 4.21 255.0 41.97326 42.51711
## 38 2002 365 365 0 4.98 289.0 49.30463 50.89087
## 39 2003 365 365 0 2.68 312.0 66.03233 67.18302
## 40 2004 366 366 0 2.77 213.0 43.52918 43.62734
## 41 2005 365 365 0 4.96 321.0 48.51677 53.08182
## 42 2006 365 365 0 2.58 333.0 70.73003 85.25151
## 43 2007 365 365 0 5.64 381.0 71.45077 67.47710
## 44 2008 366 366 0 5.30 166.0 40.58847 30.46670
## 45 2009 365 365 0 5.48 407.0 53.16282 68.14245
## 46 2010 365 365 0 4.39 365.0 64.36416 61.11077
## 47 2011 365 365 0 5.68 222.0 57.97247 47.16704
## 48 2012 366 366 0 3.55 229.0 58.51451 48.18720
# visuallize the min, max, and mean
plotdata <- melt(flow.sum,
id.var='Year',
measure.var=c("min.Q","max.Q","mean.Q","sd.Q"),
variable.name='Statistic',
value.name='Value')
ggplot(data=plotdata, aes(x=Year, y=Value))+
ggtitle("Summary statistics about Q over time")+
theme(plot.title = element_text(hjust = 0.5))+
geom_line()+
facet_wrap(~Statistic, ncol=2, scales="free_y")
# Some preliminary plots to check for outliers etc
ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=Q))+
ggtitle("Q by date")+
theme(plot.title = element_text(hjust = 0.5))+
geom_line(aes(group=Year))
ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=log(Q)))+
ggtitle("log(Q) by date")+
theme(plot.title = element_text(hjust = 0.5))+
geom_line(aes(group=Year))
ggplot(data=flow[flow$Year >= start.year,], aes(x=Date, y=log(Q)))+
ggtitle("Q by date")+
theme(plot.title = element_text(hjust = 0.5))+
geom_line()+
facet_wrap(~Year, scales="free_x")
## Warning: Removed 1461 rows containing missing values (geom_path).
# Compute the statistics on an annual basis
na.rm=list(na.rm.global=TRUE)
# Create the directory for the results of the analysis
# If you are comparing to the spreadsheet, the spreadsheet must also be in this directory.
# Similarly, if you are comparing to the HEC-SSP results, the vfa.rpt file must also be in this directory
report.dir <- file.path("08HA011","MinFlows") # directory for output files and save objects
dir.create(report.dir)
## Warning in dir.create(report.dir): '08HA011/MinFlows' already exists
cat("Reports and saved files will be found in ", report.dir, "\n")
## Reports and saved files will be found in 08HA011/MinFlows
stat.annual <- compute.Q.stat.annual(Station.code=Station.Code,
Station.Area=Station.Area,
flow=flow,
start.year=start.year,
end.year=end.year,
write.stat=TRUE, # write out statistics
write.stat.trans=TRUE, # write out statistics in transposed format
report.dir=report.dir,
na.rm=na.rm)
names(stat.annual)
## [1] "Q.stat.annual" "Q.stat.annual.trans" "dates.missing.flows" "file.stat" "file.stat.trans" "na.rm" "Date"
head(stat.annual$Q.stat.annual)
## Year SW_1_day_MIN SW_1_day_MINDOY SW_3_day_MIN SW_3_day_MINDOY SW_7_day_MIN SW_7_day_MINDOY SW_30_day_MIN SW_30_day_MINDOY SW_ANNUAL_MIN SW_ANNUAL_MAX SW_ANNUAL_MEAN SW_ANNUAL_TOTALQ
## 1 1965 3.54 276 3.746667 276 4.224286 276 5.070000 223 3.54 196 47.72625 1505094910
## 2 1966 5.80 215 5.936667 224 6.305714 225 6.445667 244 5.80 362 60.32545 1902423451
## 3 1967 4.11 243 4.110000 256 4.320000 269 4.652333 270 4.11 213 61.70992 1946083966
## 4 1968 5.10 219 5.223333 220 5.258571 223 6.569000 237 5.10 450 69.52503 2198548223
## 5 1969 4.98 228 4.980000 238 5.071429 240 5.357667 256 4.98 150 47.50655 1498166495
## 6 1970 4.25 203 4.343333 204 4.642857 206 5.035667 230 4.25 153 36.84767 1162028161
## SW_ANNUAL_YIELDMM SW_JFM_TOTALQ SW_AMJ_TOTALQ SW_JAS_TOTALQ SW_OND_TOTALQ SW_JFM_YIELDMM SW_AMJ_YIELDMM SW_JAS_YIELDMM SW_OND_YIELDMM JAN_MIN_DAILY_SW FEB_MIN_DAILY_SW MAR_MIN_DAILY_SW
## 1 1823.397 565669441 232522271 48689856 658213342 684.8298 281.5040 58.94656 796.8685 28.9 96.3 31.7
## 2 2304.754 732991673 314576353 74649600 780205824 887.3991 380.8430 90.37482 944.5591 57.2 54.4 46.7
## 3 2357.648 874480320 258498431 44503776 768601439 1058.6929 312.9521 53.87866 930.5102 106.0 58.6 49.6
## 4 2656.226 1148152321 227285567 93315456 729794880 1390.0149 275.1641 112.97271 883.5289 67.7 90.0 74.8
## 5 1815.003 443759040 540264386 76129631 438013438 537.2385 654.0731 92.16662 530.2826 30.0 28.3 36.5
## 6 1407.777 492134399 222619105 46146240 401128417 595.8044 269.5147 55.86712 485.6276 39.1 48.1 39.9
## APR_MIN_DAILY_SW MAY_MIN_DAILY_SW JUN_MIN_DAILY_SW JUL_MIN_DAILY_SW AUG_MIN_DAILY_SW SEP_MIN_DAILY_SW OCT_MIN_DAILY_SW NOV_MIN_DAILY_SW DEC_MIN_DAILY_SW JAN_MAX_DAILY_SW FEB_MAX_DAILY_SW
## 1 28.9 8.07 6.51 4.81 4.50 4.11 3.54 72.8 68.8 92.6 196.0
## 2 43.0 7.36 6.77 6.09 5.80 6.26 6.37 30.6 119.0 260.0 119.0
## 3 34.5 28.10 6.23 5.80 4.11 4.11 4.70 51.0 49.3 213.0 187.0
## 4 36.8 10.10 8.07 6.40 5.10 12.60 14.40 72.8 88.1 450.0 211.0
## 5 81.3 58.00 6.51 6.09 4.98 5.24 24.40 23.5 34.0 98.5 65.7
## 6 24.7 12.50 5.47 4.25 4.70 5.24 6.23 19.5 45.0 150.0 80.7
## MAR_MAX_DAILY_SW APR_MAX_DAILY_SW MAY_MAX_DAILY_SW JUN_MAX_DAILY_SW JUL_MAX_DAILY_SW AUG_MAX_DAILY_SW SEP_MAX_DAILY_SW OCT_MAX_DAILY_SW NOV_MAX_DAILY_SW DEC_MAX_DAILY_SW JAN_MEAN_SW FEB_MEAN_SW
## 1 86.7 70.8 56.1 15.9 9.23 11.40 9.66 67.1 171.0 188 45.10323 124.08571
## 2 117.0 140.0 43.9 31.7 26.90 13.30 8.75 43.3 182.0 362 123.84194 83.40357
## 3 138.0 77.6 39.1 32.6 7.39 6.74 6.51 167.0 156.0 187 143.38710 111.24643
## 4 135.0 76.5 39.6 20.8 16.30 16.00 19.30 122.0 141.0 177 204.73871 136.61724
## 5 103.0 141.0 79.0 64.8 12.70 6.09 36.20 39.1 49.6 150 60.32903 45.78214
## 6 72.2 138.0 23.8 11.8 5.66 5.32 11.60 27.2 71.4 153 71.65161 66.56071
## MAR_MEAN_SW APR_MEAN_SW MAY_MEAN_SW JUN_MEAN_SW JUL_MEAN_SW AUG_MEAN_SW SEP_MEAN_SW OCT_MEAN_SW NOV_MEAN_SW DEC_MEAN_SW JAN_P50_SW FEB_P50_SW MAR_P50_SW APR_P50_SW MAY_P50_SW JUN_P50_SW JUL_P50_SW
## 1 54.01613 41.90333 38.00129 8.536333 5.498065 6.489677 6.397333 33.44871 96.13000 119.27097 39.6 117.0 54.9 32.70 40.5 7.930 5.15
## 2 74.49355 73.04000 31.00968 16.281000 13.913871 6.799032 7.396667 16.87290 61.63667 214.77419 111.0 81.7 69.7 72.80 35.7 11.000 14.00
## 3 82.62581 49.75333 32.96452 15.912667 6.551613 5.538710 4.676333 70.17258 94.94667 124.90645 124.0 103.0 81.0 46.70 33.1 11.295 6.51
## 4 96.12903 53.06333 21.15484 12.764000 11.491290 7.942258 15.920000 52.24194 100.08000 123.38065 203.0 124.0 98.3 50.15 19.0 11.800 10.90
## 5 64.00000 104.38000 69.13548 32.615333 8.174516 5.536774 15.202666 31.73226 40.38667 92.71935 61.2 42.6 73.1 103.00 70.5 32.850 7.65
## 6 51.97097 60.76333 17.22903 7.320333 5.243871 5.068064 7.147667 13.30613 43.72000 94.14839 62.0 65.8 46.7 55.35 18.7 6.555 5.32
## AUG_P50_SW SEP_P50_SW OCT_P50_SW NOV_P50_SW DEC_P50_SW JAN_P80_SW FEB_P80_SW MAR_P80_SW APR_P80_SW MAY_P80_SW JUN_P80_SW JUL_P80_SW AUG_P80_SW SEP_P80_SW OCT_P80_SW NOV_P80_SW DEC_P80_SW JAN_P90_SW
## 1 5.24 5.915 25.70 95.10 111.0 33.7 101.40 38.2 29.70 34.3 7.360 4.96 4.96 4.960 19.40 79.90 89.2 29.7
## 2 6.54 7.335 7.82 54.10 182.0 91.2 59.70 54.4 51.70 23.4 6.928 6.48 6.03 6.698 6.94 35.76 165.0 66.8
## 3 5.38 4.450 73.10 98.55 118.0 118.0 77.62 64.0 39.10 30.3 7.220 6.29 5.10 4.296 47.90 65.54 86.7 113.0
## 4 6.48 15.350 38.20 89.45 121.0 90.9 106.00 79.0 39.32 13.0 8.964 8.78 5.47 13.460 21.20 80.24 106.0 76.2
## 5 5.66 7.180 32.80 42.20 104.0 36.8 34.94 40.5 89.68 64.6 6.880 6.65 5.10 5.320 25.20 36.44 40.2 33.7
## 6 5.18 5.540 9.23 48.85 88.1 47.0 60.60 40.5 37.22 12.9 6.190 4.98 4.70 5.240 7.31 22.72 62.3 44.7
## FEB_P90_SW MAR_P90_SW APR_P90_SW MAY_P90_SW JUN_P90_SW JUL_P90_SW AUG_P90_SW SEP_P90_SW OCT_P90_SW NOV_P90_SW DEC_P90_SW SW_Oct_to_Sept_TOTALQ SW_Oct_to_Sept_YIELDMM SW_ONDJFM_TOTALQ
## 1 98.68 34.3 29.40 21.6 6.940 4.96 4.96 4.779 6.23 76.06 82.7 1132277554 1370.796 1143909313
## 2 58.00 51.3 47.37 13.7 6.770 6.48 6.03 6.577 6.65 31.64 162.0 1780430969 2155.485 1391205016
## 3 68.08 57.2 36.74 29.7 6.850 6.09 4.53 4.182 7.56 54.23 85.0 1957688352 2370.083 1654686144
## 4 95.08 74.8 38.20 12.8 8.580 7.67 5.10 12.690 16.10 78.70 96.0 2237354782 2708.662 1916753760
## 5 30.95 40.5 85.50 59.7 6.645 6.51 4.98 5.240 25.20 28.90 37.4 1789947937 2167.007 1173553920
## 6 57.29 39.9 33.26 12.6 5.993 4.53 4.70 5.240 6.60 20.47 48.4 1198913181 1451.469 930147837
## SW_AMJJAS_TOTALQ SW_ONDJFM_YIELDMM SW_AMJJAS_YIELDMM
## 1 281212127 1384.878 340.4505
## 2 389225953 1684.268 471.2179
## 3 303002207 2003.252 366.8308
## 4 320601022 2320.525 388.1368
## 5 616394017 1420.767 746.2397
## 6 268765345 1126.087 325.3818
tail(stat.annual$Q.stat.annual)
## Year SW_1_day_MIN SW_1_day_MINDOY SW_3_day_MIN SW_3_day_MINDOY SW_7_day_MIN SW_7_day_MINDOY SW_30_day_MIN SW_30_day_MINDOY SW_ANNUAL_MIN SW_ANNUAL_MAX SW_ANNUAL_MEAN SW_ANNUAL_TOTALQ
## 43 2007 5.64 197 5.680000 198 5.810000 245 5.952667 259 5.64 381 71.45077 2253271393
## 44 2008 5.30 218 5.330000 220 5.448571 221 5.669000 232 5.30 166 40.58847 1283504835
## 45 2009 5.48 242 5.533333 264 5.578571 268 5.625000 271 5.48 407 53.16282 1676542751
## 46 2010 4.39 242 4.460000 260 4.561429 242 4.792000 260 4.39 365 64.36416 2029788287
## 47 2011 5.68 256 5.840000 258 5.892857 260 6.210333 264 5.68 222 57.97247 1828219678
## 48 2012 3.55 281 3.596667 282 4.138572 277 4.687667 283 3.55 229 58.51451 1850369184
## SW_ANNUAL_YIELDMM SW_JFM_TOTALQ SW_AMJ_TOTALQ SW_JAS_TOTALQ SW_OND_TOTALQ SW_JFM_YIELDMM SW_AMJ_YIELDMM SW_JAS_YIELDMM SW_OND_YIELDMM JAN_MIN_DAILY_SW FEB_MIN_DAILY_SW MAR_MIN_DAILY_SW
## 43 2729.800 1045543680 310594176 61282656 835850881 1265.7914 376.0220 74.19208 1011.9260 95.1 69.4 67.6
## 44 1550.696 555906242 292714561 77844672 357039360 673.0100 354.3760 94.24294 432.2510 53.5 42.8 46.8
## 45 2031.103 423973439 267626592 47885472 937057248 513.2850 324.0031 57.97273 1134.4519 32.8 27.3 28.5
## 46 2459.054 865581119 414793440 43291584 706122143 1047.9190 502.1712 52.41112 854.8694 99.6 64.7 53.3
## 47 2214.857 871732798 420056929 64668672 471761278 1055.3666 508.5435 78.29137 571.1396 91.7 65.1 58.1
## 48 2235.566 809023680 393163200 63979200 584203104 979.4476 475.9845 77.45666 707.2677 76.8 73.7 64.1
## APR_MIN_DAILY_SW MAY_MIN_DAILY_SW JUN_MIN_DAILY_SW JUL_MIN_DAILY_SW AUG_MIN_DAILY_SW SEP_MIN_DAILY_SW OCT_MIN_DAILY_SW NOV_MIN_DAILY_SW DEC_MIN_DAILY_SW JAN_MAX_DAILY_SW FEB_MAX_DAILY_SW
## 43 55.0 25.0 6.53 5.64 5.65 5.77 21.50 50.8 71.8 304 107.0
## 44 28.1 19.5 15.00 5.59 5.30 5.80 16.70 27.5 30.8 166 65.5
## 45 38.6 16.6 7.03 5.92 5.48 5.50 5.50 63.9 83.2 121 49.2
## 46 70.0 17.6 7.82 5.14 4.39 4.39 4.41 54.2 62.3 365 113.0
## 47 50.9 29.9 6.38 6.26 6.09 5.68 12.90 35.5 37.9 211 130.0
## 48 68.9 22.2 15.10 6.44 5.53 4.08 3.55 57.1 63.4 229 133.0
## MAR_MAX_DAILY_SW APR_MAX_DAILY_SW MAY_MAX_DAILY_SW JUN_MAX_DAILY_SW JUL_MAX_DAILY_SW AUG_MAX_DAILY_SW SEP_MAX_DAILY_SW OCT_MAX_DAILY_SW NOV_MAX_DAILY_SW DEC_MAX_DAILY_SW JAN_MEAN_SW FEB_MEAN_SW
## 43 278.0 105.0 67.4 24.5 9.42 9.06 19.20 112.0 157 381 176.93871 83.93214
## 44 71.0 47.2 83.4 49.9 15.70 14.20 18.20 32.4 111 56 101.85484 52.84483
## 45 90.2 68.0 63.3 16.4 7.03 6.50 5.76 84.2 407 239 73.64194 36.96071
## 46 99.3 133.0 66.7 67.4 7.89 5.40 9.38 85.8 108 261 176.43871 87.05357
## 47 222.0 107.0 74.5 50.4 11.50 10.40 26.80 51.7 193 135 135.17097 86.58214
## 48 132.0 94.1 85.3 25.4 17.60 7.12 5.99 68.6 130 194 128.37419 94.08966
## MAR_MEAN_SW APR_MEAN_SW MAY_MEAN_SW JUN_MEAN_SW JUL_MEAN_SW AUG_MEAN_SW SEP_MEAN_SW OCT_MEAN_SW NOV_MEAN_SW DEC_MEAN_SW JAN_P50_SW FEB_P50_SW MAR_P50_SW APR_P50_SW MAY_P50_SW JUN_P50_SW JUL_P50_SW
## 43 137.61290 70.78333 33.91935 13.99467 6.841290 7.333226 8.996000 66.63226 94.61000 153.88065 149 83.40 146.0 70.10 29.9 14.50 6.12
## 44 56.26129 36.96333 50.99355 23.27333 9.100645 6.504194 13.907667 24.57742 70.35333 40.64194 104 51.20 54.8 37.85 56.2 21.10 7.74
## 45 51.26774 54.76333 36.80645 10.45433 6.449355 5.977419 5.633333 20.92806 208.24000 127.40645 76 35.80 41.3 57.10 33.2 9.27 6.39
## 46 68.10323 85.71667 43.28065 29.58833 5.936774 4.926452 5.476667 38.01000 75.67000 152.39677 163 84.55 64.5 78.95 43.5 20.50 5.65
## 47 112.09355 77.20667 53.54839 29.51900 8.702258 7.310323 8.403000 38.41935 66.76000 73.10968 128 80.50 100.0 81.10 55.7 30.90 8.92
## 48 85.66129 79.50000 50.99355 19.49000 12.847097 6.103548 5.101000 15.77774 89.49333 115.73226 122 90.40 85.1 74.85 48.4 18.00 13.90
## AUG_P50_SW SEP_P50_SW OCT_P50_SW NOV_P50_SW DEC_P50_SW JAN_P80_SW FEB_P80_SW MAR_P80_SW APR_P80_SW MAY_P80_SW JUN_P80_SW JUL_P80_SW AUG_P80_SW SEP_P80_SW OCT_P80_SW NOV_P80_SW DEC_P80_SW
## 43 7.11 6.350 65.9 87.60 137.0 117.0 72.54 78.8 63.06 25.9 7.088 5.90 5.84 6.000 46.70 59.70 102.0
## 44 5.76 15.350 25.7 71.15 39.2 71.4 46.12 49.9 31.92 21.4 17.760 5.76 5.49 10.316 21.40 54.08 33.0
## 45 5.82 5.625 12.2 192.00 112.0 50.9 30.06 31.1 43.92 21.2 7.704 6.12 5.67 5.598 7.30 102.96 94.5
## 46 5.00 5.160 33.3 72.95 156.0 115.0 75.06 59.3 73.48 25.5 8.630 5.41 4.56 4.788 6.65 60.88 92.7
## 47 6.80 6.160 41.0 41.85 61.1 108.0 71.40 68.7 53.40 37.5 13.100 6.66 6.40 5.954 31.50 37.70 44.5
## 48 5.90 5.245 11.2 79.30 97.4 96.9 84.04 72.0 71.94 36.4 16.240 10.30 5.72 4.658 4.26 65.20 78.8
## JAN_P90_SW FEB_P90_SW MAR_P90_SW APR_P90_SW MAY_P90_SW JUN_P90_SW JUL_P90_SW AUG_P90_SW SEP_P90_SW OCT_P90_SW NOV_P90_SW DEC_P90_SW SW_Oct_to_Sept_TOTALQ SW_Oct_to_Sept_YIELDMM SW_ONDJFM_TOTALQ
## 43 110.0 72.09 72.8 56.50 25.4 6.643 5.78 5.81 5.894 22.40 54.04 89.4 2401276320 2907.114 2029399488
## 44 58.8 45.54 48.5 29.67 20.7 15.840 5.69 5.36 6.635 20.50 34.40 31.9 1762316356 2133.555 1391757123
## 45 35.0 29.03 29.7 40.16 17.0 7.522 6.11 5.62 5.578 7.17 71.59 88.3 1096524863 1327.512 781012799
## 46 111.0 70.35 57.0 73.04 22.2 8.040 5.34 4.49 4.655 4.76 57.16 70.4 2260723392 2736.953 1802638367
## 47 103.0 67.00 63.6 52.03 34.3 10.061 6.42 6.27 5.895 31.40 37.23 41.0 2062580543 2497.071 1577854941
## 48 87.4 81.48 68.2 70.95 27.6 15.890 7.43 5.61 4.148 3.98 62.52 73.1 1737927359 2104.028 1280784958
## SW_AMJJAS_TOTALQ SW_ONDJFM_YIELDMM SW_AMJJAS_YIELDMM
## 43 371876832 2456.9001 450.2141
## 44 370559234 1684.9360 448.6189
## 45 315512064 945.5361 381.9759
## 46 458085025 2182.3709 554.5824
## 47 484725601 1910.2360 586.8349
## 48 457142400 1550.5871 553.4412
head(stat.annual$Q.stat.annual.trans)
## Y1965 Y1966 Y1967 Y1968 Y1969 Y1970 Y1971 Y1972 Y1973 Y1974 Y1975 Y1976 Y1977 Y1978 Y1979 Y1980 Y1981
## SW_1_day_MIN 3.540000 5.800000 4.11 5.100000 4.980000 4.250000 6.000000 5.830000 3.11 3.280000 6.430000 4.620000 4.160000 4.64 3.840000 4.420000 4.660000
## SW_1_day_MINDOY 276.000000 215.000000 243.00 219.000000 228.000000 203.000000 223.000000 304.000000 285.00 308.000000 206.000000 223.000000 227.000000 212.00 187.000000 244.000000 252.000000
## SW_3_day_MIN 3.746667 5.936667 4.11 5.223333 4.980000 4.343333 6.293333 6.000000 4.47 3.443333 6.713333 4.976667 4.643333 4.70 3.920000 4.430000 4.833333
## SW_3_day_MINDOY 276.000000 224.000000 256.00 220.000000 238.000000 204.000000 224.000000 304.000000 253.00 309.000000 207.000000 223.000000 227.000000 213.00 188.000000 244.000000 249.000000
## SW_7_day_MIN 4.224286 6.305714 4.32 5.258571 5.071429 4.642857 6.718571 6.141429 4.57 3.712857 6.922857 5.405714 4.975714 4.78 3.941429 4.897143 4.931429
## SW_7_day_MINDOY 276.000000 225.000000 269.00 223.000000 240.000000 206.000000 225.000000 262.000000 279.00 309.000000 210.000000 225.000000 227.000000 215.00 177.000000 245.000000 253.000000
## Y1982 Y1983 Y1984 Y1985 Y1986 Y1987 Y1988 Y1989 Y1990 Y1991 Y1992 Y1993 Y1994 Y1995 Y1996 Y1997 Y1998 Y1999
## SW_1_day_MIN 5.320000 5.710000 4.120000 3.360000 3.45 3.600000 4.140000 4.080000 4.890000 4.610000 3.740000 4.89 4.800000 4.670000 4.63 7.45 3.960000 6.010000
## SW_1_day_MINDOY 213.000000 225.000000 203.000000 241.000000 297.00 221.000000 244.000000 280.000000 227.000000 233.000000 221.000000 276.00 233.000000 175.000000 239.00 231.00 264.000000 260.000000
## SW_3_day_MIN 5.360000 5.746667 4.420000 3.616667 3.49 3.733333 4.416667 4.120000 4.966667 4.706667 3.790000 5.06 4.880000 4.800000 4.65 7.52 3.996667 6.120000
## SW_3_day_MINDOY 234.000000 233.000000 203.000000 234.000000 297.00 222.000000 245.000000 281.000000 227.000000 210.000000 221.000000 247.00 234.000000 176.000000 239.00 236.00 249.000000 261.000000
## SW_7_day_MIN 5.381429 5.857143 5.021429 3.717143 3.68 3.968571 4.530000 4.368571 5.245714 4.765714 4.002857 5.09 5.037143 4.872857 4.72 7.59 4.018571 6.235714
## SW_7_day_MINDOY 219.000000 233.000000 207.000000 243.000000 298.00 314.000000 224.000000 284.000000 227.000000 214.000000 224.000000 251.00 245.000000 181.000000 243.00 236.00 250.000000 261.000000
## Y2000 Y2001 Y2002 Y2003 Y2004 Y2005 Y2006 Y2007 Y2008 Y2009 Y2010 Y2011 Y2012
## SW_1_day_MIN 5.230000 4.210000 4.980000 2.68 2.770000 4.960000 2.580000 5.64 5.300000 5.480000 4.390000 5.680000 3.550000
## SW_1_day_MINDOY 246.000000 230.000000 309.000000 271.00 228.000000 221.000000 256.000000 197.00 218.000000 242.000000 242.000000 256.000000 281.000000
## SW_3_day_MIN 5.300000 4.226667 5.106667 2.75 2.823333 5.020000 2.616667 5.68 5.330000 5.533333 4.460000 5.840000 3.596667
## SW_3_day_MINDOY 246.000000 232.000000 309.000000 277.00 229.000000 221.000000 257.000000 198.00 220.000000 264.000000 260.000000 258.000000 282.000000
## SW_7_day_MIN 5.391429 4.264286 5.317143 2.77 2.888571 5.131429 2.660000 5.81 5.448571 5.578571 4.561429 5.892857 4.138572
## SW_7_day_MINDOY 247.000000 232.000000 310.000000 277.00 233.000000 222.000000 259.000000 245.00 221.000000 268.000000 242.000000 260.000000 277.000000
head(stat.annual$dates.missing.flows)
## character(0)
stat.annual$file.stat
## [1] "08HA011/MinFlows/08HA011-annual-summary-stat.csv"
stat.annual$file.stat.trans
## [1] "08HA011/MinFlows/08HA011-annual-summary-stat-trans.csv"
# Compare the annual statistics with those in the Excel spreasheet
compare.annual <- compare.annual.stat(Q.filename=stat.annual$file.stat,
E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
report.dir=report.dir,
save.plots=TRUE,
save.comparison=TRUE)
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
names(compare.annual)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat" "plot.list" "stat.not.plotted" "file.comparison" "file.plots" "Date"
compare.annual$stats.in.Q.not.in.E
## character(0)
compare.annual$stats.in.E.not.in.Q
## [1] "SW_TOTAL_Q_CMS" "Date_50P_Cummul_Q" "Date_25P_Cummul_Q" "Date_75P_Cummul_Q"
head(compare.annual$diff.stat)
## Year Value.Q Value.E diff mean pdiff stat
## 1 1965 3.54 3.54 1.469700e-10 3.54 4.151695e-11 SW_1_day_MIN
## 2 1966 5.80 5.80 2.651399e-10 5.80 4.571378e-11 SW_1_day_MIN
## 3 1967 4.11 4.11 4.856000e-10 4.11 1.181508e-10 SW_1_day_MIN
## 4 1968 5.10 5.10 3.674296e-10 5.10 7.204503e-11 SW_1_day_MIN
## 5 1969 4.98 4.98 -7.348966e-11 4.98 1.475696e-11 SW_1_day_MIN
## 6 1970 4.25 4.25 0.000000e+00 4.25 0.000000e+00 SW_1_day_MIN
names(compare.annual$plot.list)
## [1] "plot.daily" "plot.yieldmm" "plot.mean" "plot.per" "plot.wyear"
l_ply(compare.annual$plot.list, function(x){plot(x)})
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
# Compute the long-term statistics
stat.longterm <- compute.Q.stat.longterm(Station.code=Station.Code,
Station.Area=Station.Area,
flow=flow,
start.year=start.year,
end.year=end.year,
write.stat=TRUE, # write out statistics
write.stat.trans=TRUE,
report.dir=report.dir) # write out statistics in transposed format
names(stat.longterm)
## [1] "Q.stat.longterm" "Q.stat.longterm.trans" "file.stat" "file.stat.trans" "na.rm" "Date"
head(stat.longterm$Q.stat.longterm)
## Month mean median max min
## 1 1 110.46828 98.65 450.0 15.60
## 2 2 92.69174 82.90 303.0 14.70
## 3 3 77.38508 68.85 359.0 15.80
## 4 4 57.23731 52.40 154.0 8.67
## 5 5 35.27665 31.00 121.0 6.47
## 6 6 18.18226 13.60 73.4 2.99
tail(stat.longterm$Q.stat.longterm)
## Month mean median max min
## 8 8 6.069792 5.550 87.9 2.77
## 9 9 8.989847 6.055 84.9 2.58
## 10 10 26.823743 14.200 312.0 2.72
## 11 11 85.113521 71.100 407.0 3.28
## 12 12 113.266129 104.000 425.0 24.20
## 13 13 53.262478 37.400 450.0 2.58
head(stat.longterm$Q.stat.longterm.trans)
## JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC Longterm
## mean 110.4683 92.69174 77.38508 57.23731 35.27665 18.18226 9.739046 6.069792 8.989847 26.82374 85.11352 113.2661 53.26248
## median 98.6500 82.90000 68.85000 52.40000 31.00000 13.60000 6.755000 5.550000 6.055000 14.20000 71.10000 104.0000 37.40000
## max 450.0000 303.00000 359.00000 154.00000 121.00000 73.40000 64.800003 87.900002 84.900002 312.00000 407.00000 425.0000 450.00000
## min 15.6000 14.70000 15.80000 8.67000 6.47000 2.99000 2.850000 2.770000 2.580000 2.72000 3.28000 24.2000 2.58000
stat.longterm$file.stat
## [1] "08HA011/MinFlows/08HA011-longterm-summary-stat.csv"
stat.longterm$file.stat.trans
## [1] "08HA011/MinFlows/08HA011-longterm-summary-stat-trans.csv"
# Compare the longterm statistics with those in the Excel spreasheet
compare.longterm <- compare.longterm.stat(Q.filename=stat.longterm$file.stat,
E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
report.dir=report.dir,
save.comparison=TRUE,
save.plots=TRUE)
## Warning: Removed 2 rows containing missing values (geom_point).
names(compare.longterm)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat" "plot.list" "stat.not.plotted" "file.comparison" "file.plots" "Date"
compare.longterm$stats.in.Q.not.in.E
## character(0)
compare.longterm$stats.in.E.not.in.Q
## character(0)
head(compare.longterm$diff.stat)
## Month Value.Q Value.E diff mean pdiff stat
## 1 1 110.46828 110.46830 -2.042926e-05 110.46829 1.849333e-07 mean
## 2 2 92.69174 92.69174 3.771534e-07 92.69174 4.068900e-09 mean
## 3 3 77.38508 77.38508 6.253091e-07 77.38508 8.080487e-09 mean
## 4 4 57.23731 57.23731 -4.460229e-06 57.23731 7.792520e-08 mean
## 5 5 35.27665 35.27665 -3.460968e-06 35.27665 9.810931e-08 mean
## 6 6 18.18226 18.18226 3.881854e-06 18.18226 2.134968e-07 mean
names(compare.longterm$plot.list)
## [1] "plot.allstat"
l_ply(compare.longterm$plot.list, function(x){plot(x)})
## Warning: Removed 2 rows containing missing values (geom_point).
# Compute the long-term percentile statistics
percentile.longterm <- compute.Q.percentile.longterm(Station.code=Station.Code,
Station.Area=Station.Area,
flow=flow,
start.year=start.year,
end.year=end.year,
write.stat=TRUE, # write out statistics
write.stat.trans=TRUE, # write out statistics in transposed format
report.dir=report.dir)
names(percentile.longterm)
## [1] "Q.percentile.stat" "Q.percentile.stat.trans" "file.stat" "file.stat.trans" "na.rm" "Date"
head(percentile.longterm$Q.percentile.stat)
## Month P01 P02 P05 P10 P15 P20 P25 P30 P35 P40 P45 P50 P55 P60 P65 P70 P75 P80 P85 P90 P95 P99
## 1 Jan 321.000 282.26 239.650 194.00 169.900 152.60 141.000 129.00 120.550 112.00 106.000 98.65 92.215 84.94 78.2000 70.81 65.800 58.600 49.60 40.41 31.480 21.8740
## 2 Feb 245.450 229.80 186.000 154.00 133.000 122.00 114.000 107.00 101.000 95.40 89.150 82.90 77.275 73.10 68.5250 64.10 60.000 54.900 49.60 43.90 36.950 23.2550
## 3 Mar 209.520 188.26 161.650 132.30 112.950 101.00 92.900 85.40 79.000 74.84 71.385 68.85 65.115 62.60 58.9450 55.40 51.300 47.200 41.91 36.07 29.235 20.3960
## 4 Apr 131.610 120.00 107.000 93.91 85.000 77.34 72.425 68.00 63.700 59.70 56.000 52.40 48.655 45.66 43.6000 41.30 38.575 35.800 32.30 28.48 23.490 13.4000
## 5 May 84.082 79.13 72.800 65.83 60.295 54.30 49.000 43.70 40.055 36.62 33.400 31.00 29.000 26.40 23.7000 21.31 18.350 16.100 14.20 12.20 9.867 7.4122
## 6 Jun 63.954 57.20 49.005 40.40 32.215 26.82 23.025 20.63 18.800 16.84 15.200 13.60 12.000 10.70 9.0965 8.16 7.580 6.948 6.55 6.28 5.580 3.9500
tail(percentile.longterm$Q.percentile.stat)
## Month P01 P02 P05 P10 P15 P20 P25 P30 P35 P40 P45 P50 P55 P60 P65 P70 P75 P80 P85 P90 P95 P99
## 7 Jul 40.500 35.504 25.3650 17.500 14.395 12.300 10.80 9.799 8.5700 7.754 7.107 6.755 6.493 6.250 6.0790 5.850 5.6600 5.464 5.2500 4.960 4.3535 3.0900
## 8 Aug 13.313 11.700 9.6495 8.096 7.179 6.716 6.43 6.120 5.9400 5.800 5.670 5.550 5.460 5.368 5.2545 5.140 5.0000 4.860 4.6100 4.350 3.8800 3.0274
## 9 Sep 37.705 34.476 23.3150 16.700 14.200 11.920 9.35 8.263 7.3705 6.718 6.270 6.055 5.860 5.660 5.5000 5.340 5.2300 5.030 4.7900 4.390 3.8895 2.8639
## 10 Oct 157.130 125.000 90.9650 64.090 49.790 39.520 32.50 27.200 22.4000 19.000 15.985 14.200 12.800 11.580 10.2000 8.765 7.6175 6.674 5.7805 5.267 4.4235 3.1861
## 11 Nov 304.000 282.220 213.0500 161.100 134.150 120.000 111.00 103.000 95.5000 87.500 79.600 71.100 64.800 57.200 51.8950 46.900 42.5000 37.700 31.4000 20.090 10.7950 4.7238
## 12 Dec 306.520 269.520 231.0000 195.000 173.000 157.600 146.00 135.000 126.0000 118.000 111.000 104.000 94.900 88.100 81.3000 73.710 66.0750 58.600 51.4050 44.270 36.5000 26.7740
head(percentile.longterm$Q.percentile.stat.trans)
## JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC
## P01 321.00 245.45 209.52 131.61 84.082 63.954 40.500 13.3130 37.705 157.130 304.00 306.52
## P02 282.26 229.80 188.26 120.00 79.130 57.200 35.504 11.7000 34.476 125.000 282.22 269.52
## P05 239.65 186.00 161.65 107.00 72.800 49.005 25.365 9.6495 23.315 90.965 213.05 231.00
## P10 194.00 154.00 132.30 93.91 65.830 40.400 17.500 8.0960 16.700 64.090 161.10 195.00
## P15 169.90 133.00 112.95 85.00 60.295 32.215 14.395 7.1790 14.200 49.790 134.15 173.00
## P20 152.60 122.00 101.00 77.34 54.300 26.820 12.300 6.7160 11.920 39.520 120.00 157.60
percentile.longterm$file.stat
## [1] "08HA011/MinFlows/08HA011-longterm-percentile-stat.csv"
percentile.longterm$file.stat.trans
## [1] "08HA011/MinFlows/08HA011-longterm-percentile-stat-trans.csv"
# Compare the longterm percentile statistics with those in the Excel spreasheet
compare.percentile.longterm <- compare.percentile.longterm.stat(Q.filename=percentile.longterm$file.stat,
E.filename=file.path(report.dir,"08HA011_STREAMFLOW_SUMMARY.xlsx"),
save.comparison=TRUE,
save.plots=TRUE,
report.dir=report.dir)
names(compare.percentile.longterm)
## [1] "stats.in.Q.not.in.E" "stats.in.E.not.in.Q" "diff.stat" "plot.list" "stat.not.plotted" "file.cmparsion" "file.plots" "Date"
compare.percentile.longterm$stats.in.Q.not.in.E
## character(0)
compare.percentile.longterm$stats.in.E.not.in.Q
## character(0)
head(compare.percentile.longterm$diff.stat)
## Month Value.Q Value.E diff mean pdiff stat
## 1 Apr 131.610 131.610 1.136868e-13 131.610 8.638161e-16 P01
## 2 Aug 13.313 13.313 -4.482903e-10 13.313 3.367312e-11 P01
## 3 Dec 306.520 306.520 4.547474e-13 306.520 1.483581e-15 P01
## 4 Feb 245.450 245.450 -5.684342e-14 245.450 2.315886e-16 P01
## 5 Jan 321.000 321.000 0.000000e+00 321.000 0.000000e+00 P01
## 6 Jul 40.500 40.500 0.000000e+00 40.500 0.000000e+00 P01
names(compare.percentile.longterm$plot.list)
## [1] "plot.allstat"
l_ply(compare.percentile.longterm$plot.list, function(x){plot(x)})
# Compute the volume frequency analysis (similar to HEC SSP)
vfa.analysis <- compute.volume.frequency.analysis( Station.Code="08HA011", flow,
start.year=1965, end.year=2016, use.water.year=FALSE,
roll.avg.days=c(1,3,7,15),
use.log=FALSE,
use.max=FALSE,
fit.distr="PIII",
write.stat=TRUE,
write.plotdata=TRUE,
write.quantiles=TRUE,
report.dir=report.dir)
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
## Warning in min(x$Q, na.rm = TRUE): no non-missing arguments to min; returning Inf
names(vfa.analysis)
## [1] "start.year" "end.year" "use.water.year" "use.max" "roll.avg.days" "Q.stat" "Q.stat.trans"
## [8] "plotdata" "prob.plot.position" "freqplot" "fit.distr" "fit" "fitted.quantiles" "fitted.quantiles.trans"
## [15] "file.stat" "file.stat.trans" "file.plotdata" "file.quantile" "file.quantile.trans" "date"
vfa.analysis$file.stat
## [1] "08HA011/MinFlows/08HA011-annual-vfa-stat.csv"
head(vfa.analysis$Q.stat)
## Year Measure value
## 1 1965 Q001-avg 3.540000
## 2 1965 Q003-avg 3.746667
## 3 1965 Q007-avg 4.224286
## 4 1965 Q015-avg 4.633333
## 5 1966 Q001-avg 5.800000
## 6 1966 Q003-avg 5.936667
head(vfa.analysis$Q.stat.trans)
## Year Q001-avg Q003-avg Q007-avg Q015-avg
## 1 1965 3.54 3.746667 4.224286 4.633333
## 2 1966 5.80 5.936667 6.305714 6.419333
## 3 1967 4.11 4.110000 4.320000 4.415333
## 4 1968 5.10 5.223333 5.258571 5.563333
## 5 1969 4.98 4.980000 5.071429 5.262000
## 6 1970 4.25 4.343333 4.642857 4.863333
vfa.analysis$freqplot
ggsave(plot=vfa.analysis$freqplot,
file=file.path(report.dir, paste(Station.Code,"-annual-vfa-frequency-plot.png",sep="")), h=6, w=6, units="in", dpi=300)
vfa.analysis$fitted.quantiles.trans
## distr prob Q001-avg Q003-avg Q007-avg Q015-avg
## 1 PIII 0.010 2.416806 2.581846 2.675951 2.838811
## 2 PIII 0.050 2.988016 3.145622 3.285016 3.450660
## 3 PIII 0.100 3.308550 3.462336 3.620229 3.787091
## 4 PIII 0.200 3.711903 3.861207 4.036056 4.204141
## 5 PIII 0.500 4.530930 4.672124 4.862202 5.031834
## 6 PIII 0.800 5.413735 5.547480 5.729219 5.899323
## 7 PIII 0.900 5.901163 6.031295 6.198934 6.368843
## 8 PIII 0.950 6.317356 6.444656 6.595474 6.764989
## 9 PIII 0.975 6.688462 6.813419 6.945789 7.114789
## 10 PIII 0.980 6.800741 6.925022 7.051201 7.220017
## 11 PIII 0.990 7.131889 7.254258 7.360602 7.528800
HEC.filename <- file.path(report.dir,"HEC-vfa-min.rpt")
compare.with.HEC <- compare.frequency.with.hec(Q.file.stat=vfa.analysis$file.stat,
Q.file.plotdata=vfa.analysis$file.plotdata,
Q.file.quantile=vfa.analysis$file.quantile,
HEC.filename=HEC.filename,
report.dir=report.dir,
save.comparison=TRUE,
save.plots=TRUE)
names(compare.with.HEC)
## [1] "stats.in.Q.not.in.HEC" "stats.in.HEC.not.in.Q" "diff.stat" "plot.list" "stat.not.plotted" "file.comparison" "file.plots" "Date"
compare.with.HEC$stats.in.Q.not.in.HEC
## character(0)
compare.with.HEC$stats.in.HEC.not.in.Q
## [1] "Q030-avg" "Q060-avg" "Q090-avg" "Q120-avg" "Q183-avg" "Q030-pp" "Q060-pp" "Q090-pp" "Q120-pp" "Q183-pp" "Q030-avg-q" "Q060-avg-q" "Q090-avg-q" "Q120-avg-q" "Q183-avg-q"
head(compare.with.HEC$diff.stat)
## Year Measure value.Q value.HEC diff mean pdiff
## 1 0.01 Q001-avg-q 2.416806 2.438 -0.02119401 2.427403 0.008731144
## 2 0.01 Q003-avg-q 2.581846 2.544 0.03784636 2.562923 0.014766873
## 3 0.01 Q007-avg-q 2.675951 2.612 0.06395081 2.643975 0.024187369
## 4 0.01 Q015-avg-q 2.838811 2.738 0.10081146 2.788406 0.036153798
## 5 0.05 Q001-avg-q 2.988016 2.976 0.01201559 2.982008 0.004029361
## 6 0.05 Q003-avg-q 3.145622 3.108 0.03762155 3.126811 0.012031925
names(compare.with.HEC$plot.list)
## [1] "plot.avg" "plot.pp" "plot.quant"
l_ply(compare.with.HEC$plot.list, function(x){plot(x)})
# show the differences between the MLE and moment estimates (?) from HEC
compare.with.HEC$diff.stat[ grepl('-q$', compare.with.HEC$diff.stat$Measure),]
## Year Measure value.Q value.HEC diff mean pdiff
## 1 0.01 Q001-avg-q 2.416806 2.438 -0.0211940053 2.427403 0.0087311441
## 2 0.01 Q003-avg-q 2.581846 2.544 0.0378463612 2.562923 0.0147668730
## 3 0.01 Q007-avg-q 2.675951 2.612 0.0639508081 2.643975 0.0241873688
## 4 0.01 Q015-avg-q 2.838811 2.738 0.1008114576 2.788406 0.0361537980
## 5 0.05 Q001-avg-q 2.988016 2.976 0.0120155866 2.982008 0.0040293612
## 6 0.05 Q003-avg-q 3.145622 3.108 0.0376215530 3.126811 0.0120319251
## 7 0.05 Q007-avg-q 3.285016 3.228 0.0570163190 3.256508 0.0175084220
## 8 0.05 Q015-avg-q 3.450660 3.378 0.0726602881 3.414330 0.0212809790
## 9 0.10 Q001-avg-q 3.308550 3.290 0.0185499803 3.299275 0.0056224414
## 10 0.10 Q003-avg-q 3.462336 3.433 0.0293360484 3.447668 0.0085089539
## 11 0.10 Q007-avg-q 3.620229 3.579 0.0412285130 3.599614 0.0114535920
## 12 0.10 Q015-avg-q 3.787091 3.740 0.0470913506 3.763546 0.0125124961
## 13 0.20 Q001-avg-q 3.711903 3.693 0.0189030938 3.702452 0.0051055614
## 14 0.20 Q003-avg-q 3.861207 3.847 0.0142065088 3.854103 0.0036860737
## 15 0.20 Q007-avg-q 4.036056 4.020 0.0160556821 4.028028 0.0039859909
## 16 0.20 Q015-avg-q 4.204141 4.192 0.0121414722 4.198071 0.0028921552
## 17 0.50 Q001-avg-q 4.530930 4.530 0.0009301536 4.530465 0.0002053108
## 18 0.50 Q003-avg-q 4.672124 4.689 -0.0168756631 4.680562 0.0036054778
## 19 0.50 Q007-avg-q 4.862202 4.893 -0.0307980917 4.877601 0.0063141885
## 20 0.50 Q015-avg-q 5.031834 5.074 -0.0421657228 5.052917 0.0083448277
## 21 0.80 Q001-avg-q 5.413735 5.437 -0.0232648779 5.425368 0.0042881662
## 22 0.80 Q003-avg-q 5.547480 5.576 -0.0285196961 5.561740 0.0051278369
## 23 0.80 Q007-avg-q 5.729219 5.766 -0.0367809254 5.747610 0.0063993431
## 24 0.80 Q015-avg-q 5.899323 5.938 -0.0386773799 5.918661 0.0065348189
## 25 0.90 Q001-avg-q 5.901163 5.932 -0.0308372016 5.916581 0.0052119965
## 26 0.90 Q003-avg-q 6.031295 6.047 -0.0157054324 6.039147 0.0026006043
## 27 0.90 Q007-avg-q 6.198934 6.207 -0.0080657874 6.202967 0.0013003112
## 28 0.90 Q015-avg-q 6.368843 6.367 0.0018431866 6.367922 0.0002894487
## 29 0.95 Q001-avg-q 6.317356 6.348 -0.0306438498 6.332678 0.0048390033
## 30 0.95 Q003-avg-q 6.444656 6.435 0.0096555098 6.439828 0.0014993429
## 31 0.95 Q007-avg-q 6.595474 6.559 0.0364736611 6.577237 0.0055454383
## 32 0.95 Q015-avg-q 6.764989 6.704 0.0609888210 6.734494 0.0090561841
## 33 0.99 Q001-avg-q 7.131889 7.140 -0.0081113549 7.135944 0.0011366898
## 34 0.99 Q003-avg-q 7.254258 7.154 0.1002581369 7.204129 0.0139167602
## 35 0.99 Q007-avg-q 7.360602 7.174 0.1866024313 7.267301 0.0256769915
## 36 0.99 Q015-avg-q 7.528800 7.282 0.2467995701 7.405400 0.0333269745